/*//////////////////////////////////////////////////////////////////////////////

Copyright Laurent Excoffier, Genetics and Biometry Lab, University of Geneva
E-mail: laurent.excoffier@anthro.unige.ch
URL: http://anthropologie.unige.ch/~laurent

The use of the following source code in only allowed if it is kept intact
with its original banner

Extensively modified by Christian Anderson for the Hadly lab, Stanford
Now at Scripps Institution of Oceanography, senka@ucsd.edu

//////////////////////////////////////////////////////////////////////////////*/

#include "genealgy.h"
#include "public.h"
#include "time.h"
#include "migrmat.h"
#include "deme.h"
#include "cstring.h"
#include <fstream.h>
#include "coalmain.h"
#include "cond_var.h"

typedef MY_TArrayAsVector<my_string> NameArray;
float UNIT_TIME;

struct TInput {
	 int num_pop;
    int num_samp;
    int num_loci;
    Mut_Type data_type;
    double mut_rate;
    double new_mut_rate; // Added by Rumen Kostadinov - global mutation rate change 
    long time_of_mut_rate_change; // Added by Rumen Kostadinov - the mut. rate change occurrs at this time (in generations)
    double transition_rate;
    double a_c; // Added by Rumen Kostadinov 12 parameter evo model
    double a_g; // Added by Rumen Kostadinov
    double a_t; // Added by Rumen Kostadinov
    double c_a; // Added by Rumen Kostadinov
    double c_g; // Added by Rumen Kostadinov
    double c_t; // Added by Rumen Kostadinov
    double g_a; // Added by Rumen Kostadinov
    double g_c; // Added by Rumen Kostadinov
    double g_t; // Added by Rumen Kostadinov
    double t_a; // Added by Rumen Kostadinov
    double t_c; // Added by Rumen Kostadinov
    double t_g; // Added by Rumen Kostadinov
  double alpha; // Kostadinov 8-26-07 LOH event prob. P(0->1) 
  double beta; // Kostadinov 8-26-07 reversion P(1->0)
    double gamma_a;
    double mut_ratio;
    double prop_sites;
    double num_rate_categories;
    int range_constraint;
};

//Internal functions
bool LoadFile(char *path, TDemeCollection *myDemes, TEventArray *myEvents,
       TDemeSizeArray *myDemeSizes, TGrowthArray *myGrowthRates, TSampArray *mySamples,
       TMigrMatArray *myMigrationMatrices, TInput *t, TNode *first_node);
void GetFiles(int argc, char* argv[], NameArray *infile_name, int *num_rand_samples, int *num_files);

int main(int argc, char *argv[]) {
   bool trace=1;

   NameArray infile_name(1,5);
   int num_rand_samples=1, num_files=0, i;

   printf("+----------------------------+\n|       SERIAL SIMCOAL       |\n| by the Hadly Lab, Stanford |\n|       via SimCoal 1.0      |\n+----------------------------+");
   
   //Handle the input
   GetFiles(argc, argv, &infile_name, &num_rand_samples, &num_files);
   printf("\n# of files=%d",num_files);
   //Reads and generate genealogies for the required number of files
   for (int nf=0; nf<num_files; ++nf) {
      printf("\nNow reading %s",infile_name[nf].c_str());
   	//Creating input and output file names
		my_string outfile_name, arl_batch_file_name,
	             arl_generic_name, paup_batch_file_name,
	             paup_generic_name, Stat_file_name,
	             paup_name, log_file_name,
	             paup_true_trees_name, paup_mut_trees_name, temp_in;

		outfile_name		= infile_name[nf]+".gen";
		arl_batch_file_name = infile_name[nf]+".arb";
		arl_generic_name	= infile_name[nf];
		paup_generic_name	= infile_name[nf];
		paup_batch_file_name= infile_name[nf]+".bat";
		paup_name			= infile_name[nf]+".paup";
		paup_true_trees_name= infile_name[nf]+"_true_trees.trees";
		paup_mut_trees_name	= infile_name[nf]+"_mut_trees.trees";
		log_file_name 		= infile_name[nf]+".log";
		Stat_file_name		= infile_name[nf]+"_stat.csv";
		temp_in				= infile_name[nf]+".par";
		
		//Opening result files
		ofstream ofs(outfile_name.c_str());
		
		#ifdef _PRINT_ARLEQUIN_OUTPUT_
			ofstream ArlBatch(arl_batch_file_name.c_str());
		#endif

      #ifdef _PRINT_SEQ_PAUP_BATCH_
        ofstream PaupBatch(paup_batch_file_name.c_str());
      #endif

      cout << "\nOutput file name : " << outfile_name << "\n";
      if (!ofs) {
   	  cout << "\nUnable to open the freaking output file (" << outfile_name << ")\n";
        return 1;
	   }

   	//Initialization of a bunch of object and vectors
   TDemeCollection   myDemes;
 	TEventArray       myEvents(10,10);  
 	TDemeSizeArray    myDemeSizes(10,10);
 	TGrowthArray      myGrowthRates(10,10);
	TSampArray 	      mySamples(10,10); //Christian 6/29/04
   TMigrMatArray     myMigrationMatrices(1,10);  //Loro_16_9_99
   TNode				   first_node;
   TInput t;
    
    //Christian: Read all the input
    if(!LoadFile(temp_in.c_str(),&myDemes,&myEvents,&myDemeSizes,&myGrowthRates,
        &mySamples,&myMigrationMatrices,&t,&first_node)) return 1;
    if(trace) cout<<"\nInitializing Demes";
   	//Initialization of the demes
   	myDemes.create_demes(&myMigrationMatrices);
      myDemes.initialize_deme_size(myDemeSizes);
   	myDemes.initialize_growth_rates(myGrowthRates);
   	myDemes.initialize_events(myEvents);
   	myDemes.create_lineages(mySamples); //Christian 06/29/04
   if(trace) cout<<"\nWriting stuff to outfile.";
    //Output initial conditions
   	ofs 	<< "\nDemes initial conditions"
         	<< "\n========================\n"
         	<< myDemes
         	<< endl;
   	//Initialize the generator of random numbers once for all
   	time_t start=time(NULL);
   	long lidum=-start; ran3(&lidum);
   	cout << "\nRandom generator initialized with : " << start << endl;
      ofs  << "Random generator initialized with : " << start << "\n\n";

   	cout << "\nBuilding " << num_rand_samples << " genealogies ...\n";

   #ifdef _PRINT_SEQ_PAUP_
      ofstream PaupFile(paup_name.c_str());
   #endif

   //Christian: stat file
   ofstream StatInfile(Stat_file_name.c_str());
   if(!StatInfile) { cout<<"Could not open statfile!"; continue; }
   if(t.data_type==DNA) {
     	for(i=0;i<=myDemes.stat_grps();i++) {
    		StatInfile << "GROUP "<<i<<",Haptypes,PrivHaps,SegSites,PairDiffs,HapDiver,NucltdDiv,TajimasD,MismatDist,";
            for(int j=i+1; j<=myDemes.stat_grps();j++)
 		      StatInfile <<i<< " VS "<<j<<",PrivTo"<<i<<",PrivTo"<<j<<",PairDiffs,MeanDiv(Hs-bar),PoolDiv(Ht),Fst,";
        }
        StatInfile << "COMBINED,Haptypes,PrivHaps,SegSites,PairDiffs,HapDiver,NucltdDiv,TajimasD,MismatDist,MRCA\n";
    } else {
        for(i=0;i<=myDemes.stat_grps();i++) {
            StatInfile << "SAMPLE "<<i+1<<",Alleles,PrivAlleles,Heterozyg,AllelicVar,";
            for(int j=i+1; j<=myDemes.stat_grps();j++)
                StatInfile <<i<< " VS "<<j<<",PrivTo"<<i<<",PrivTo"<<j;
        }
        StatInfile << "COMBINED,Alleles,PrivAlleles,Heterozyg,AllelicVar\n";
    }

   #ifdef _PRINT_TREE_PAUP_
   	ofstream
         PaupTrueTreesFile(paup_true_trees_name.c_str()),
         PaupMutTreesFile(paup_mut_trees_name.c_str());
   #endif
    int tot_num_nodes;
    for (i=0; i<t.num_samp; ++i)
        tot_num_nodes+=mySamples[i].n();

   #ifdef _PRINT_SEQ_PAUP_
	   PaupFile
         << "#NEXUS\n"
      	 << "\n[Simulated data generated by program sim_coal.exe (Laurent Excoffier)"
         << "\n modified for ancient DNA by Liz Hadly's Lab, Stanford Univ."
         << "\n\tSimulation conditions:"
         << "\n\tNo. of taxa   : " << tot_num_nodes
         << "\n\tNo. of loci   : " << t.num_loci
         << "\n\tNo. of demes  : " << t.num_pop
         << "\n\tMutation rate : " << (t.mut_rate/t.num_loci) << " per locus (site) per generation"
         << "\n\tAlpha         : " << t.gamma_a << " (0=homogeneity)"
         << "\n\t%Transitions  : " << (100*t.transition_rate) << " (ts/tv=";
      if (t.transition_rate!=1)
      	PaupFile << (t.transition_rate/(1-t.transition_rate));
      else
			PaupFile << "inf";
      PaupFile << ")\n]\n\n";  
      write_PAUP_header(PaupFile,log_file_name.c_str());
  	#endif

   #ifdef _PRINT_TREE_PAUP_
      write_PAUP_trees_section_header(PaupTrueTreesFile);
      write_PAUP_trees_section_header(PaupMutTreesFile);
  	#endif


      float mut_rate_per_locus=0.0;
      if (t.num_loci) mut_rate_per_locus=t.data_type/t.num_loci;
      Mut_Model mut_mod;

      if (fabs(t.gamma_a)>1e-7) mut_mod=K80_GAMMA;
      else mut_mod=K80_NOGAMMA;
      
      tree_output_type tree_type;   

    //Build unequal mutation rates if needed
    //Create the vector of gamma distributed mutation rates once for all

    //Loro_26_2_99
    if (!t.gamma_a==0.0 && t.data_type==DNA) {
     	cout << "\nInitializing unequal mutation rates...";
   	first_node.mut_rates.get_new_rates(t.num_loci,t.gamma_a, t.mut_ratio, t.prop_sites, (int)t.num_rate_categories);
     	cout << "done\n\n";
    }
	
    my_float tree_mut_length=0.0, tree_mut_length2=0.0, tree_exp_length=0.0, tree_exp_length2=0.0;
    int tot_cases=0;
    for (; tot_cases<num_rand_samples; ++tot_cases) {
   	   if(trace) cout<<"\nRandom Sample " << tot_cases <<"\n - Building tree...";
   		if (myDemes.build_tree()) { //HERE WE GO!!
         	//Add mutations
         	if(trace) cout<<"done.\n - Sprinkling Mutations...";

		// Modified by Rumen Kostadinov
		/* original:
      		long mut=myDemes.sprinkle_mutations(t.mut_rate, t.num_loci, t.data_type,
                    t.gamma_a, t.mut_ratio, t.prop_sites, t.transition_rate, t.range_constraint);
		*/
		long mut=myDemes.sprinkle_mutations
		  (t.mut_rate, t.new_mut_rate, t.time_of_mut_rate_change, t.num_loci, t.data_type,
		   t.gamma_a, t.mut_ratio, t.prop_sites, t.transition_rate, 
		   t.a_c, t.a_g, t.a_t, t.c_a, t.c_g, t.c_t, 
		   t.g_a, t.g_c, t.g_t, t.t_a, t.t_c, t.t_g, t.alpha, t.beta, t.range_constraint);
		
      		if(trace) cout<<"done.";
            tree_mut_length+=mut;
            tree_mut_length2+=mut*mut;
            tree_exp_length+=first_node.tree_exp_length();
            tree_exp_length2+=first_node.tree_exp_length()*first_node.tree_exp_length();

         #ifdef _PRINT_SEQ_PAUP_
			//Writing data to PAUP file
			if(trace) cout<<"\n - Writing to PAUP file...";
            write_PAUP_replicate_number((tot_cases+1), PaupFile);
            write_PAUP_data_header(tot_num_nodes, t.num_loci, t.data_type, PaupFile);
            myDemes.write_samples_to_PAUP_file(PaupFile, t.data_type);
            write_PAUP_end_matrix(PaupFile);
			   write_PAUP_tree_header(PaupFile, (tot_cases+1));
            tree_type=MUT_RATE;
            myDemes.print_gene_tree(PaupFile, tree_type, mut_rate_per_locus);
            write_PAUP_end(PaupFile);
            if (t.num_loci) {
            	PaupFile
              		<< "\n\t[Tree length = " << (first_node.tree_exp_length()/t.num_loci)
                    << "\n\tNumber of sites hit by mutations = "  << first_node.count_polym_sites()
                    << "]\n";
            }
            write_PAUP_block(PaupFile,mut_mod);
			if(trace) cout<<"done.";
         #endif

         #ifdef _PRINT_TREE_PAUP_
				if(trace) cout<<"\n - Writing to PAUP Tree...";
				write_PAUP_tree_name(PaupTrueTreesFile, (tot_cases+1));
         	write_PAUP_tree_name(PaupMutTreesFile, (tot_cases+1));
            tree_type=GENERATIONS;
            myDemes.print_gene_tree(PaupTrueTreesFile, tree_type, mut_rate_per_locus);
            tree_type=NUM_MUT;
            myDemes.print_gene_tree(PaupMutTreesFile, tree_type, mut_rate_per_locus);
            if(trace) cout<<"done.";
         #endif

            if (!fmod((tot_cases+1),10))
         		cout << "\nGenealogy # " << (tot_cases+1) << "/" << num_rand_samples;

         
   		#ifdef _PRINT_SEQ_PAUP_BATCH_
			if(trace) cout<<"\n - Writing to PAUP Sequence...";
			//Build paup file name               
            char PaupOut[150], PaupLog[150];
            sprintf(PaupOut,"%s_%d.pau",paup_generic_name.c_str(),tot_cases);
            sprintf(PaupLog,"%s_%d.log",paup_generic_name.c_str(),tot_cases);
            //Open paup output file
      		ofstream PaupOutFile(PaupOut);
            
            //Write name in paup batch file
            PaupBatch << "paup " << PaupOut << "\n";
            PaupOutFile
         		<< "#NEXUS\n"
      			<< "\n[Simulated data generated by program sim_coal.exe (Laurent Excoffier)"
         		<< "\n\tSimulation conditions:"
         		<< "\n\tNo. of taxa   : " << tot_num_nodes
         		<< "\n\tNo. of loci   : " << t.num_loci
         		<< "\n\tNo. of demes  : " << t.num_pop
         		<< "\n\tMutation rate : " << (t.data_type/t.num_loci) << " per locus (site) per generation"
         		<< "\n\tAlpha         : " << t.gamma_a << " (0=homogeneity)"
         		<< "\n\t%Transitions  : " << (100*t.transition_rate) << " (ts/tv=";
            if (t.transition_rate!=1) PaupOutFile << (t.transition_rate/(1-t.transition_rate));
      		else	PaupOutFile << "inf";
      		PaupOutFile << ")\n]\n\n";

            write_PAUP_header(PaupOutFile,PaupLog);
            write_PAUP_replicate_number((tot_cases+1), PaupOutFile);
            write_PAUP_data_header(tot_num_nodes, t.num_loci, t.data_type,PaupOutFile);
            myDemes.write_samples_to_PAUP_file(PaupOutFile, t.data_type);
				write_PAUP_end_matrix(PaupOutFile);
				write_PAUP_tree_header(PaupOutFile, (tot_cases+1));
            tree_type=MUT_RATE;
            myDemes.print_gene_tree(PaupOutFile, tree_type, mut_rate_per_locus);
            write_PAUP_end(PaupOutFile);
            	PaupOutFile
            		<< "\n\t[Tree length = " << (first_node.tree_exp_length()/t.num_loci)
                  << "\n\tNumber of sites hit by mutations = "  << first_node.count_polym_sites()
                  << "]\n";
				write_PAUP_block(PaupOutFile,mut_mod);
				write_PAUP_footer(PaupOutFile);  
            if(trace) cout<<"done.";
      	#endif

        #ifdef _PRINT_ARLEQUIN_OUTPUT_
            //Build Arlequin file name
			if(trace) cout<<"Writing ARLEQUIN file...";
            char ArlOut[150];
            sprintf(ArlOut,"%s_%d.arp",paup_generic_name.c_str(),tot_cases);

            //Open Arlequin output file
      		ofstream ArlInFile(ArlOut);

            //Write name in Arlequin batch file
            ArlBatch << ArlOut << "\n";

            ArlInFile
            	<< "#Arlequin input file written by program the simulation program sim_coal.exe\n\n"
      			<< "#Simulation parameters:\n"
               << "#======================\n"
               << "\n#Deme sizes\n";

            for (i=0; i<t.num_pop; ++i) {
   				ArlInFile << "#Deme #" << i << "\t" << myDemeSizes[i] << "\n";
   			}
            ArlInFile << "\n#Sample sizes\n";
   			for (i=0; i<t.num_samp; ++i) {
   				ArlInFile << "#Deme #" << mySamples[i].deme() << "\t" << mySamples[i].n() << "\n";
   			}
            ArlInFile << "\n#Growth rates\n";
   			for (i=0; i<t.num_pop; ++i) {
   				ArlInFile << "#Deme #" << i << "\t" << myGrowthRates[i] << "\n";
   			}
            for (i=0; i<myMigrationMatrices.GetItemsInContainer(); ++i) {
            	ArlInFile
            		<< "\n#Migration matrix " << i << "\n"
   					<< myMigrationMatrices[i];
            }
            ArlInFile << "\n\n#Historical events\n";
   			for (i=0; i<myEvents.GetItemsInContainer(); ++i) {
   				ArlInFile << "#Event #" << i << "\n" << myEvents[i] << "\n";
   			}            
      		ArlInFile
		  << "\n#Data type							 : ";
		switch (t.data_type) {
		case DNA     	: ArlInFile << "DNA (Transition rate="
					    << setprecision(3)
					    << (t.transition_rate*100.0)
					    << "%)\n";
		  break;
		case LOH        : ArlInFile << "LOH (alpha="
					    << setprecision(3)
					    << (t.alpha)
					    << " beta= "
					    << setprecision(3)
					    << (t.beta)
					    << ")\n";
		  break;
		case RFLP      : ArlInFile << "RFLP\n"; break;
		case MICROSAT  : {
		  if (t.range_constraint)
		    ArlInFile
		      << "MICROSAT"
		      << " (Range constraint : " << t.range_constraint
		      << " [" << (10000+first_node.min_mic) << ";"
		      << (10000+first_node.max_mic) << "])\n";
		  else       
		    ArlInFile
		      << "MICROSAT"
		      << " (No range constraint)\n";
		  
		}; break;
		}
		ArlInFile
		  //            	<< setprecision(10)
            	<< "\n#Mutation rate per generation : " << t.data_type
                << "\n#Number of loci to simulate   : " << t.num_loci
             	<< "\n#Gamma parameter              : " << t.gamma_a
//            	<< "\n#Mutation ratio               : " << t.mut_ratio
               << "\n\n";

            //Loro_20_9_99: Get the number of demes with non-empty samples
            int num_real_samples=0;
            for (int d=0; d<t.num_pop; ++d) {
            	if (myDemes[d].sample_size()) ++num_real_samples;
            }

            write_Arlequin_header(num_real_samples, t.data_type, ArlInFile);
			   myDemes.write_samples_to_Arlequin_file(ArlInFile, t.data_type);
            myDemes.write_group_section_to_Arlequin_file(ArlInFile);
            ArlInFile.close();
			if(trace) cout<<"done.";
    	#endif
         if(trace) cout<<"\n - Writing to Statfile...";
       	//Christian: calculate and run statistics
      	myDemes.write_Statfile(StatInfile, t.data_type);
    	   if(trace) cout<<"done.";
      	} 
      	myDemes.reset(myDemeSizes,myGrowthRates);
      	if(trace) system("PAUSE");
      } //end random coal sample
      if (t.num_loci) {
         tree_exp_length/=t.num_loci;
         tree_exp_length2/=t.num_loci*t.num_loci;
      }
      else {
         tree_exp_length=0;
         tree_exp_length2=0;
      }

    #ifdef _PRINT_SEQ_PAUP_
        write_PAUP_footer(PaupFile);   
    #endif
   
   #ifdef _PRINT_TREE_PAUP_
      write_PAUP_end(PaupTrueTreesFile);
      write_PAUP_end(PaupMutTreesFile);
      PaupTrueTreesFile
        << "\nMean length of the trees : "
        << (tree_exp_length/num_rand_samples);
      if (num_rand_samples>1)
      	PaupTrueTreesFile
            << " +- "
         	<< setiosflags(ios::showpoint | ios::fixed | ios::right)
         	<< setw(10)
         	<< setprecision(5)
            << ( sqrt((tree_exp_length2 - tree_exp_length*tree_exp_length/num_rand_samples)/
            			 (num_rand_samples-1)) )
            << endl;

      PaupMutTreesFile
        << "\nMean length of the trees : "
        << (tree_mut_length/num_rand_samples);
      if (num_rand_samples>1)
      	PaupMutTreesFile
            << " +- "
         	<< setiosflags(ios::showpoint | ios::fixed | ios::right)
         	<< setw(10)
         	<< setprecision(5)
            << ( sqrt((tree_mut_length2 - tree_mut_length*tree_mut_length/num_rand_samples)/
            			 (num_rand_samples-1)) )
            << endl;
   #endif
        
   	start=time(NULL)-start;
   	cout << "\n\ndone in " <<  start << "s\n" << endl;
	
   	//Print to output file

      ofs   << "\nFinal results"
   			<< "\n=============\n";

      ofs 	<< "\nData type : ";
      switch (t.data_type) {
      	case DNA 		: ofs << "DNA"; break;
      	case RFLP 		: ofs << "RFLP"; break;
      	case MICROSAT 	: ofs << "MICROSAT"; break;
      	case LOH 		: ofs << "LOH"; break;
      }

      if (t.data_type==DNA) {
      	if (t.gamma_a>0.0) {
         	ofs	<< "\nHeterogeneity of mutation rates: Gamma value = "
            	<< setiosflags(ios::showpoint | ios::fixed | ios::right)
         		<< setw(10)
         		<< setprecision(5)
         		<< t.gamma_a;
         }
         else ofs << "\nNo heterogeneity of mutation rates";
         ofs    << "\nProportion of transitions : "
            	<< setiosflags(ios::showpoint | ios::fixed | ios::right)
         		<< setw(10)
         		<< setprecision(5) << (100*t.transition_rate) << " (ts/tv=";
         if (t.transition_rate!=1) ofs << (t.transition_rate/(1-t.transition_rate)) << ")";
      	 else ofs << "inf";
      }
            
      ofs   << "\nNumber of loci : " << t.num_loci;
      ofs	<< "\nTotal mutation rate per generation : " << t.data_type;
      ofs	<< "\nMutation rate per locus per generation : " << t.data_type/t.num_loci;
	
      ofs
      	<< "\nMean length of the trees (subst. rate per site) : " 
         	<< setiosflags(ios::showpoint | ios::fixed | ios::right)
         << (tree_exp_length/num_rand_samples);
      if (num_rand_samples>1)
      	ofs << " +- "
         	<< setiosflags(ios::showpoint | ios::fixed | ios::right)
         	<< setw(10)
         	<< setprecision(5)
            << ( sqrt((tree_exp_length2 - tree_exp_length*tree_exp_length/num_rand_samples)/
            			 (num_rand_samples-1)) )
            << endl;

      ofs
      	<< "\nMean length of the trees (total number of mutations) : "
         << (tree_mut_length/num_rand_samples);
      if (num_rand_samples>1)
      	ofs << " +- "
         	<< setiosflags(ios::showpoint | ios::fixed | ios::right)
         	<< setw(10)
         	<< setprecision(5)
            << ( sqrt((tree_mut_length2 - tree_mut_length*tree_mut_length/num_rand_samples)/
            			 (num_rand_samples-1)) )
            << endl;
      ofs   << "\n\nTotal time for " << tot_cases
         	<< " simulations: " << start << " ms"
         	<< endl;

      //Do not forget to reset the static member of TNode
      TNode resetNode;
      resetNode.reset_node_count();
	}
   if(trace) system("PAUSE");
   return 1;  
}

//------------------------------------------------------------------------------
int
TMeanMat::compute_mean(const int& n) {
	_num_updates=n;
	if (!_num_updates) return 0;
	for (int i=0; i<_size; ++i) {
   	for (int j=0; j<_size; ++j) {
      	_data[i][j]/=_num_updates;
      }
   }
   return 1;
};

//Compute the s.d. from the values present in the matrix, assuming that
//those are the sum of square values
int
TMeanMat::compute_sd(const int& n, const TMeanMat& mean) {
   if (n<2) return 0;
	_num_updates=n   ;
	if (!_num_updates) return 0;
	for (int i=0; i<_size; ++i) {
   	for (int j=0; j<_size; ++j) {
      	double s2=_data[i][j],
         		 m=mean._data[i][j]; //assumes that mean contains mean values
      	_data[i][j]=sqrt( (s2-n*m*m)/(n-1) );
      }
   }
   return 1;
};

int
TMeanMat::update_with(const TMigrationMatrix& MM) {
	if (!MM.size() || _size!=MM.size()) return 0;
   for (int i=0; i<_size; ++i) {
   	for (int j=0; j<_size; ++j) {
      	_data[i][j]+=MM(i,j);
      }
   }return 1;
};

int
TMeanMat::update_with_square(const TMigrationMatrix& MM) {
	if (!MM.size() || _size!=MM.size()) return 0;
   for (int i=0; i<_size; ++i) {
   	for (int j=0; j<_size; ++j) {
      	_data[i][j]+=MM(i,j)*MM(i,j);
      }
   }return 1;
};

bool LoadFile(char *path, TDemeCollection *myDemes, TEventArray *myEvents,
       TDemeSizeArray *myDemeSizes, TGrowthArray *myGrowthRates, TSampArray *mySamples,
       TMigrMatArray *myMigrationMatrices, TInput *t, TNode *first_node) {
   //Initialize input file
   ifstream ifs;
   ifs.open(path);
  	if (!ifs) {
      cout << "Unable to open input file (2)(" << path << ")" << endl;
      return 0;
   }
   bool is_anc_file=0;
   my_string line(20);

   line.read_to_delim(ifs);   //Reads a blank line

   BayesRead(&(t->num_pop),&ifs,1);        //Reads the number of populations to simulate
   cout<<"# populations:"<<t->num_pop<<endl;
/*Christian 6/29/04: If the rest of the first line contains 'with ancient', then
  not all the samples are the same age, and appropriate action will be taken below.
  If not, the program is still backwards compatible */
   line.read_to_delim(ifs);
   if(line.contains("with ancient"))
   {   is_anc_file=1;
       cout<<"Contains ancient sampling information\n";
   } else cout<<"No ancient samples\n";
  	//Read deme sizes
  	line.read_to_delim(ifs);   //Reads a blank line
	long cur_size; int i;
  	for (i=0; i<t->num_pop; ++i) {
     	BayesRead(&cur_size,&ifs,2);
     	line.read_to_delim(ifs);
  		myDemeSizes->Add(cur_size);
  	}
  	cout << "\nDeme sizes\n";
  	for (i=0; i<t->num_pop; ++i) cout << "Deme " << i << "\t" << (*myDemeSizes)[i] << "\n";
  	cout << "\n";

   //Reads the sample sizes
  	line.read_line(ifs);  //Reads a blank line
/*Christian 6/29/04: If the file contains ancient dates, then the format here
  is 'samp age deme stat_grp'. The number of sample groups is on the first line. Otherwise
  one sample group per deme is assumed */
   if(is_anc_file) { ifs>>t->num_samp; line.read_to_delim(ifs);
   } else t->num_samp=t->num_pop; //no priors
  	for (i=0; i<t->num_samp; ++i) {
      TSamp cursamp;
      if(is_anc_file) ifs >> cursamp; //supports priors
      else {
         BayesRead(&cur_size,&ifs,2);
         cursamp.SetVals(cur_size,0,i,0);
      }
      line.read_to_delim(ifs);
   	mySamples->Add(cursamp);
	}
   cout << "\nSamples\n----------";
	for (i=0; i<t->num_samp; ++i)
		cout << "\nSample " << i << " - " << (*mySamples)[i];
	cout << endl;
//Christian 6/29/04: End modifications

   line.read_line(ifs);  //Reads a blank line
   float cur_growth;
   for (i=0; i<t->num_pop; ++i) {
      BayesRead(&cur_growth,&ifs,3);
      line.read_to_delim(ifs);
      myGrowthRates->Add(cur_growth);
   }
   cout << "\nGrowth rates\n";
   for (i=0; i<t->num_pop; ++i) cout << "Deme" << i << ":\t" << (*myGrowthRates)[i] << "\n";
   cout << endl;

   //Loro_16_9_99
   //Reading how many migration matrices to read
   int numMigMat=0; //Default
   line.read_line(ifs);  //Reads a blank line
   ifs >> numMigMat; //no priors
   line.read_line(ifs);

   //Loro_16_9_99
   if (numMigMat) {
      for (int m=0; m<numMigMat; ++m) {
         //Reading migration matrix
         TMigrationMatrix myMigrationRates(t->num_pop);
   		line.read_line(ifs);  //Reads a blank line
   		for (i=0; i<t->num_pop; ++i)
   			for (int j=0; j<t->num_pop; ++j)
      			BayesRead(&(myMigrationRates(i,j)),&ifs,3);
         myMigrationMatrices->Add(myMigrationRates);
   		cout << "\nMigration matrix\n";
   		cout << (*myMigrationMatrices)[m];
      	line.read_to_delim(ifs); //Reads until the end of line
   	}
   } else { //Then a model without migration is assumed
   	TMigrationMatrix myMigrationRates(t->num_pop); //Enough to create an ^matrix filled with zeroes
      myMigrationMatrices->Add(myMigrationRates);
   }

   	int num_events;
   	//Read historical events
  	   line.read_line(ifs); //Reads a blank line
   	ifs >> num_events;   //Reads the number of historical events to read, no priors
   	line.read_line(ifs); //Reads the string " historical events\n"
   	for (i=0; i<num_events; ++i) {
   		THistoricalEvent curevent;
   		ifs >> curevent;  //rewritten to support priors
   		(*myEvents).Add(curevent);
   	}
   	cout << "\nHistorical events\n";
   	for (i=0; i<num_events; ++i) {
   		cout << "Event " << i << "\n" << (*myEvents)[i] << "\n";
   	}
   	cout << endl;

      //Loro_26_2_99 to avoid a bug when no historical event is defined
      if (num_events) line.read_to_delim(ifs); //Reads until the end of line

   	  //Read mutation rate
      line.read_line(ifs); //Reads a blank line
      
      // Modified by Rumen Kostadinov
      // original: ifs >> t->mut_rate;
      ifs >> t->mut_rate >> t->new_mut_rate >> t->time_of_mut_rate_change; // incorporate rate change at time t
      cout << "\nMutation rate per generation : " << t->mut_rate << 
	" Mut. rate transitions to: " << t->new_mut_rate << " at time: " << t->time_of_mut_rate_change;
      line.read_to_delim(ifs); //Reads until the end of line

      //Read number of loci
      line.read_line(ifs);
      ifs >> t->num_loci;
      cout << "\nNumber of loci to simulate : " << t->num_loci;
      line.read_to_delim(ifs); //Reads until the end of line

      //Reads data type
      t->transition_rate=0.0;
      // Added by Rumen Kostadinov
      t->a_c = 0.0; t->a_g = 0.0; t->a_t = 0.0;
      t->c_a = 0.0; t->c_g = 0.0; t->c_t = 0.0;
      t->g_a = 0.0; t->g_c = 0.0; t->g_t = 0.0;
      t->t_a = 0.0; t->t_c = 0.0; t->t_g = 0.0;
      // Kostadinov 8-26-07
      t->alpha = 0.0; t->beta = 0.0;
      //Loro_26_2_99
      t->range_constraint=0; //A value of zero assumes no range constraint
      line.read_line(ifs);  //Read comment line
      bool par=line.read_to_delim(ifs, ' ');
      line.to_upper();
      //Loro_26_2_99!
      if (line.contains("MICROSAT")) {
      	t->data_type=MICROSAT;
      	if(par) ifs >> t->range_constraint; //Extract range constraint for microsat data
      	line.read_to_delim(ifs);
      } else if (line.contains("DNA")) {
      	t->data_type=DNA;
	// Modified by Rumen Kostadinov
	// original:
        // if(par) ifs >> t->transition_rate; //Extract transition rate for DNA
	if(par) ifs >> t->transition_rate
		    >> t->a_c >> t->a_g >> t->a_t
		    >> t->c_a >> t->c_g >> t->c_t
		    >> t->g_a >> t->g_c >> t->g_t
		    >> t->t_a >> t->t_c >> t->t_g; //Extract transition rate, 12 parameters for DNA
      	line.read_to_delim(ifs);
      } else if (line.contains("LOH")) { // Kostadinov 8-26-07
	t->data_type=LOH;
        if(par) ifs >> t->alpha >> t->beta; //Extract alpha and beta probs.
        line.read_to_delim(ifs);
      } else if (line.contains("RFLP")) t->data_type=RFLP;
      // Loro_28_2_99 : Be careful: only one node allowed to be created before
      // building the tree
      first_node->reset_node_count(); //Do not forget to reset the node count...

      if (t->data_type==MICROSAT && t->range_constraint) {//Fix minimum and maximum size for microsats
         double y=2.0, x=t->range_constraint;
         if (fmod(x,y)==0.0) { //Even number
         	first_node->min_mic=-(t->range_constraint/2)+1;
            first_node->max_mic=t->range_constraint/2;
         } else {
         	first_node->min_mic=(1-t->range_constraint)/2;
            first_node->max_mic=(t->range_constraint+1)/2; //bug?!
         }
      }

      cout << "\nData type : ";
      switch (t->data_type) {
      case DNA: if(t->transition_rate>0.0000001){
	  cout << "DNA\n" << "Transition rate : " << t->transition_rate << ", Using 1 parameter DNA evo model\n";
	} else {
	  cout << "DNA\n" << "Transition rate : " << t->transition_rate << ", Using 12 parameter DNA evo model\n";
	}
	break;
      case LOH: cout << "LOH\n"; 
	break;
      case RFLP: cout << "RFLP\n"; break;
      case MICROSAT: 
         	if (t->range_constraint)
         		cout << "MICROSAT"
          			 << "\nRange constraint : " << t->range_constraint
                     << " (" << (10000+first_node->min_mic) << ";"
                     << (10000+first_node->max_mic) << ")\n";
         	else
         		cout << "MICROSAT" << "\n No range constraint\n";
            break;
         default : break;
      }

    //Loro_26_2_99
    t->gamma_a=0.0; t->mut_ratio=1.0; t->prop_sites=1.0; t->num_rate_categories=0;
    if (t->data_type==DNA) {
   	//Reads gamma parameter
   	line.read_line(ifs);
   	ifs >> t->gamma_a;
   	cout << "\nGamma parameter : " << t->gamma_a;
   	//If gamma parameter is negative, then it means that we have a
   	//two mutation rate model. So read the mutation ratio and the proportion
   	//of sites with the high mutation rate
   	if (t->gamma_a<0.0) ifs >> t->mut_ratio >> t->prop_sites;
   	else if (t->gamma_a>0.0) ifs >> t->num_rate_categories;
   }

   //loro_27_2_99
   ifs.close();
   return 1;
}

void GetFiles(int argc, char* argv[], NameArray *infile_name, int *num_rand_samples, int *num_files) {
   my_string usage;
   usage ="\nsim_coal.exe usage:\n";
   usage+="\n [1] : sim_coal <no arguments>";
   usage+="\n        The program will prompt you for an input file name to process and the";
   usage+="\n        number of simulations to perform";
   usage+="\n [2] : sim_coal <batch file name>";
   usage+="\n        The program will open the batch file that should contain, on the first";
   usage+="\n        line, the number of simulations to perform, and on the following lines";
   usage+="\n        the names of the input files to process";
   usage+="\n [3] : sim_coal <file name> <number of simulations to perform>\n";

   my_string batchFile;
   ifstream ifs_batch;

   my_string Path="";

   if (argc==2) {
   	  //Then one supposes that one receives a file containing the name of
      //the files to be processed
      batchFile=argv[1];
      Path=extract_path(batchFile);
      ifs_batch.open(batchFile.c_str());
      if (!ifs_batch) {
        cout << "Unable to open input file (" << batchFile << ")\n";
        cout << usage; return;
      } else {
     	my_string line;
      	//Number of random samples to draw must be on first line
      	ifs_batch >> *num_rand_samples;
        line.read_to_delim(ifs_batch);
      	while (ifs_batch.good()) {
     	    line.read_to_delim(ifs_batch);
            if (!line.is_null()) {
          	    ++(*num_files);
                infile_name->Add(remove_extension(line));
            }
      	}
        ifs_batch.close();
      }
   } else if (argc==3) {
   	//Then one supposes that one receives a file name to be processed
    //and the number of permutations to be performed
      *num_files=1;
      batchFile=argv[1];
      Path=extract_path(batchFile);
      ifs_batch.open(batchFile.c_str());
      if (!ifs_batch) {
   	  cout << "Unable to open input file (1) (" << batchFile << ")" << endl;
        cout << usage; return;
      }
      else {
         *num_rand_samples=atoi(argv[2]);
         ifs_batch.close();
         infile_name->Add(remove_extension(batchFile));
         if (!(*num_rand_samples)) cout << usage;
      }
   } else {
      *num_files=1;
      my_string curName;
      cout << "\n\n[Note: do not use spaces in the path name]\nGeneric input file name  (w/o .par extention)  : ";
      cin >> curName;
      infile_name->Add(curName);
      cout << "\n\nNo. of random samples to generate : ";
      cin >> *num_rand_samples;
   }
   printf("\nDone handling command paramaters");
}

